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SUMMARY 


We  derive  from  Maxwell’s  equations  a  set  of  three  coupled  first-order  differential  equations 
for  high-frequency  light  ray  trajectories  in  a  medium  whose  index  of  refraction  n(r)  at  any  point  r 
is  a  function  only  of  r.  The  derived  equations,  displayed  in  Equations  33  or  34  and  involving  the 
variables  shown  in  Figure  2,  are  simple  and  singularity-free  when  r  is  bounded  away  from  zero. 
These  equations  should  therefore  be  especially  well  suited  to  numerical  ray  tracing  in  spherically 
symmetric  models  of  Earth’s  atmosphere.  We  obtain  analytical  solutions  to  these  ray  equations 
for  two  mathematically  simple  cases,  and  we  discuss  strategies  for  obtaining  numerical  solutions 
for  more  realistic  cases. 


MAXWELL’S  EQUATIONS 


In  the  absence  of  free  electric  charges  and  currents,  Maxwell’s  equations  read 


VxE  =  -a(B, 

(la) 

V  •  D  =  0, 

(lb) 

V  •  B  =  0, 

(lc) 

VxH  =  dt  D, 

(Id) 

where  V  ■  d/dr  and  df  *  d/dt.  In  a  linear,  isotropic,  stationary  medium,  we  have  the  constitutive 
relations 


D  =  £  E,  H  =  B /fi,  (2) 

where  e  and  n  are  scalars  that  are  independent  of  t  but  may  depend  on  r.  Substituting  Equations 
2  into  Equations  lb  and  Id,  and  noting  that 


and 


V  •  (cE)  =  e  V  •  E  +  Ve  •  E, 


Vx(B//*)  =  x  B  +  x 
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we  obtain 


V  x  E  =  -  dtB, 

V  •  E  =  -e-We-E, 

V  •  B  =  0, 

V  x  B  =  ye  d(E  —  /iV(p  '')xB. 


(3a) 

(3b) 

(3c) 

(3d) 


Equations  3  define  the  E  and  B  fields  in  any  region  that  has  no  free  electric  charges  or  currents, 
and  that  is  characterized  by  scalar  point  functions  e(r)  and  /dr). 


PLANE  WAVE  SOLUTIONS  FOR  CONSTANT  e  AND  p 


In  the  special  case  in  which  e  and  y.  are  independent  of  r,  Maxwell’s  equations  3  evidently 
simplify  to 


V  x  E  =  -a,B, 

(4a) 

V-E  =  0, 

(4b) 

V  •  B  =  0, 

(4c) 

TxB  =  /it  dtE. 

(uations  admit  solutions  of  the  form 

(4d) 

E0  exp(i£0[nk  •  r  -  c<] ), 

(5a) 

Bq  exp(i£0[nk  •  r  -  ct] ), 

(5b) 

where  c  is  the  speed  of  light  in  vacuum,  and  where  the  constant  scalar  n  and  the  constant  vectors 
Eo,  Bo,  and  It  (a  circumflex  denotes  a  vector  of  unit  length)  satisfy  the  following  three  conditions: 


n  =  c(ye)m\ 

A 

{k,  E0,  Bq}  is  a  right-orthogonal  triad; 
B0  =  nE^c. 


(6a) 

(6b) 

(6c) 


Before  presenting  the  proof,  let  us  observe  that  the  functions  in  Equations  5  describe  a  plane  wave 
disturbance  that  propagates  in  the  direction  £  with  speed  c/n.  The  wave  is  plane  because  the 
surfaces  of  constant  phase  or  "wave  fronts”  are  defined  by  ft  •  r  =  constant,  which  is  the  equation 
of  a  plane  normal  to  it.  The  wave  is  evidently  sinusoidal  with  period  T  and  wavelength  A  such 
that  kQc  =  2n/T  and  kQn  =  2n/X\  thus,  the  frequency  v=  l/T  and  wavelength  A  of  the  sinusoidal 
plane  wave  described  by  Equations  5  are 


v  =  kQc/2n, 


A  =  2n/k0n. 


(7a) 

(7b) 


To  prove  that  Equations  5  satisfy  Equations  4  if  Equations  6  hold,  we  let  <p  denote  the 
argument  of  the  exponentials  in  Equations  5,  and  we  note  that  the  divergence  and  curl  of  E(r,i) 
have  the  form 


V  *  E(r,<)  =  T(e*)  *  E0  =  e*  V<p  *  E0, 

where  denotes  either  "  •  ”  or  "  x  ”  throughout.  Since  the  gradient  of  <p  is  given  by 

V<p  =  ikQn  V(ic  •  r)  =  ikQnii, 

then  we  get 

Y*E(r,fl  =  e^(iJfe0n)i c  *  E„. 

It  similarly  follows  that 

V*B(r,«  =  e*(i*0n)lc*B0. 


Substituting  these  divergence  and  curl  relations  into  Equations  4  and  then  dividing  through  by 
e*,  we  obtain 


ikQn  k  x  Eq  =  +  ikQc  BQ, 

(8a) 

ikQn  k  •  Eg  =  0, 

(8b) 

ik0n  k  •  B0  =  0, 

(8c) 

ikQn  k  x  B0  =  ye(  -  ikQc)  EQ. 

(8d) 

The  middle  two  equations  show  thatlc  must  be  perpendicular  to  both  E0  and  BQ,  and  the  first  (or 
last)  equation  then  shows  that  EQ  and  Bq  must  be  perpendicular  to  each  other,  with  1c,  EQ,  B0 
forming  a  right-orthogonal  triad,  thus  establishing  Equation  6b.  Given  Equation  6b,  the  first  and 
last  of  Equations  8  become 

nEQ  =  cB0  and  nBQ  =  yecEQ, 

respectively;  together,  these  two  equations  imply  Equations  6a  and  6c. 


DISTORTED  PLANE  WAVE  SOLUTIONS  FOR  SPATIALLY 
DEPENDENT  cANDp 


If  the  scalars  e  and  y.  that  characterize  the  electromagnetic  properties  of  the  medium  are 
functions  of  position  r,  then  we  must  use  Maxwell’s  equations  in  the  form  of  Equations  3  instead 
of  Equations  4.  Let  us  investigate  the  possibility  of  a  solution  to  Equations  3  of  the  following  form 
[cf.  Equations  5): 


E(r,<)  =  EQ(r)  exp(  ikQ[L( r)  -  c<] ), 
B(r ,t)  =  B0(r)exp(  ik0[U r)  -  cf] ). 


(9a) 

(9b) 
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These  two  functions  describe  a  distorted  plane  wave  that  propagates  locally  in  the  direction  of 
VLi r)  with  speed  c/|VL(r)|.  Notice  that  VL(r)  is  normal  to  the  (generally  nonplanar)  surface  L(r) 
=  constant,  the  surface  of  constant  phase  or  "wave  front.”  At  any  fixed  point  r  the  wave  is 
temporally  sinusoidal  with  frequency  v  satisfying  kQc  =  2 nv.  Although  the  wave  is  not  in 
general  spatially  sinusoidal,  we  can  ascribe  to  it  a  local  wavelength  .l(r)  satisfying 
fc0|VL(r)|  =  2n/A(r).  Hence,  the  distorted  plane  waves  in  Equations  9  have  frequency  and  local 
wavelength  [cf.  Equations  7] 

v  =  k^c/2n,  (10a) 

.l(r)  =  2n/fc()|VL(r)|.  (10b) 

To  investigate  whether  E(r,t)  and  B(r,t)  in  Equations  9  can  satisfy  Maxwell’s  equations  (3), 
we  again  let  cp  denote  the  phase  of  the  waves,  and  we  note  that 

V  •  E(r,f)  =  Vie*)  *  EQ  +  e*V  *  EQ 
=  e*V<t>*  E0  +  e*V  *  EQ 
=  e*likQVL)*  E0  +  e*V*E0, 

where  again  denotes  either  or  "x”  throughout.  Substituting  this  and  the  analogous 
relation  for  V  *  B(r,t)  into  Equations  3  and  then  dividing  through  by  e^,  we  get 


ikQVL  x  E0  +  V  x  E0  =  +  ikQc  Bq,  (11a) 

ikQVL  •  E0  +  V  •  E0  =  -e_1Vc  •  E0,  (lib) 

ik0VL  •  B0  +  V  •  B0  =  0,  (11c) 

ik0^L x  BQ  +  V  x  B„  =  pe(  -  tAQc)E0  -  pV(^ - 1 )  x  BQ.  (lid) 


In  general,  then,  if  EQ(r),  B0(r),  Z-(r),  and  are  chosen  to  satisfy  Equations  11,  the 
distorted  plane  waves  of  Equations  9  will  satisfy  Maxwell’s  equations.  But  now  we  consider  the 
special  case  of  geometrical  optics,  which  corresponds  to  light  waves  of  very  high  frequencies  v.  or 
by  Equations  10a,  light  waves  with  very  large  kQ.  In  that  limit,  we  can  evidently  approximate 
Equations  11  by  dropping  all  terms  that  do  not  contain  a  factor  kQ.  But  Equations  11  thus 
approximated  are  seen  to  be  identical  to  the  plane  wave  Equations  8,  provided  we  replace 

(12) 

We  have  previously  shown  that  satisfaction  of  Equations  8  requires  satisfaction  of  Equations  6 
Therefore,  using  Equation  12,  we  conclude  that  satisfaction  of  Equations  11  in  the  limit  of 
infinitely  large  k Q  requires: 


|TL(r)|  =  e  ip(r)  dr))1/2; 

(13a) 

{VL(r),E0(r),  BQ(r)}  is  a  right-orthogonal  triad: 

(13b) 

BQ(r)  =  |VL(r)|  £0(r)/c. 

(13c) 

In  summary,  if  Equations  13  are  satisfied,  then  the  distorted  plane  waves  E(r,t)  and  B(r,f) 
in  Equations  9  satisfy  Maxwell’s  equations  (3)  in  the  limit  k0-*&. 
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w 
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We  define  the  index  of  refraction  n  as  the  scalar  point  function 

n  =  n(r)  *  c  (/dr)  e(r))1/2.  1 14) 

This  definition  agrees  with  the  plane  wave  equation  (6a).  The  requirement  in  Equation  13a  now 
takes  the  form 

|VL(r)|  =  n(r).  (15) 

This  equation  is  called  the  eikonal  equation,  and  is  a  key  result  of  this  section.  As  we  have 
previously  noted,  the  distorted  plane  waves  in  Equations  9  propagate  with  speed  c/|VL(r)|,  so  we 
may  conclude  from  Equation  15  that  n(r)  is  the  ratio  of  the  propagation  speed  in  a  vacuum  to  the 
propagation  speed  in  the  medium  at  point  r.  This  interpretation  of  n(r)  is  in  agreement  with  the 
customary  definition  of  the  index  of  refraction. 

We  shall  also  define  the  unit  ray  uectork  as  the  vector  point  function 

$  ■$(*■)  *  VL(r)/|VL(r)|.  (16) 

In  close  analogy  with  the  interpretation  of  1c  in  the  plane  wave  case,1c(r)  defines  the  direction  of 
propagation  of  the  distorted  plane  waves  in  Equations  9  at  the  point  r;  equivalently,  1c(  r)  is  the 
unit  tangent  to  the  light  ray  at  r.  With  Equation  16,  the  eikonal  equation  (15)  can  be  written 

VL(r)  =  n(r)1c(r),  (17) 

thus  making  formal  the  correspondence  noted  earlier  in  Equation  12. 


EQUATION  FOR  THE  TURNING  RATE  OfIc 


We  shall  now  use  the  eikonal  equation  to  derive  an  equation  for  the  local  turning  rate  of 
the  unit  ray  vector  1c.  To  this  end,  we  consider  a  distorted  plane  wave  [cf.  Equations  91  with 
eikonal  function  L(r)  in  the  neighborhood  of  some  point  P.  Let  n(P)  and  Vn(P)  denote  the 
respective  values  of  the  index  of  refraction  and  its  gradient  at  P,  and  let  L(P)  and  VHP)  denote  the 
respective  values  of  the  eikonal  function  and  its  gradient  at  P  The  direction  of  propagation  of  the 
wave  front  at  P,  i.e.,  the  tangent  to  the  ray  trajectory  at  P,  is  ic(P)  =  VL(P)/|VL(P)|.  Let  G  be  the 
unit  vector  in  the  plane  oflc(P)  and  Vn(P)  that  makes  a  right  angle  withlc(P)  and  an  acute  angle 
with  Vn(P),  as  shown  in  Figure  1.  And  let  P'  be  a  point  whose  position  vector  relative  to  P  is  Udp, 
where  dp  is  a  positive  infinitesimal.  Since  (jt  by  construction  is  normal  to  VHP),  then  the 
smallness  of  dp  ensures  that  every  point  on  the  line  segment  PP'  has  L(r)  =  L(P).  thus,  the 
segment  PP'  lies  in  the  wave  front  through  P. 


Now  let  Q  be  a  point  whose  position  vector  relative  to  P  is  1c(P)ds,  where  ds  is  another 
positive  infinitesimal,  also  shown  in  Figure  l.  As  we  move  along  the  light  ray  the  distance  ds 
from  P  to  Q,  the  value  of  the  eikonal  function  L(r)  changes  by  an  amount 

dL  =  VL(P) -1c(P)ds  =  n(P)1c(P)-1c(P)ds  =  n(P)ds,  (18) 
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where  the  second  equality  has  invoked  the  eikonal  equation  (17).  Similarly,  as  we  move  an 
infinitesimal  distance  ds'  in  the  ray  direction  it(P')  =it(P)  from  point  P'  to  some  point  Q'  (Figure 
1),  the  eikonal  function  changes  by  an  amount 

d L'  =  n(P’)ds.  (19) 


Let  us  choose  ds'  so  that 


dL'  =  dL. 


(20) 


In  that  case,  every  point  on  the  line  segment  QQ'  has  L(r)  =  L(P)  +  dL;  thus,  the  segment  QQ' 
lies  in  the  wave  front  through  Q. 


Given  Equation  20,  we  may  deduce  from  Equations  18  and  19  that 

n(P)ds  =  n(P')ds'.  (21) 


Since  the  angle  between  {I  and  Vn(P)  is  acute,  then  n(P')  >  n(P),  so  it  follows  from  Equation  21 
that  ds' S  ds.  Evidently,  in  moving  along  the  light  ray  the  infinitesimal  distance  ds  from  P  to  Q, 
the  wave  front  and  its  attendant  unit  normal  have  rotated  through  the  infinitesimal  angle 


(ds  —  ds')  ds  /  ds'  \ 

de  =  =  57 1 1  ~  57  )• 


(22) 


the  sense  of  this  rotation  being  from  k(P)  towards  u  Using  Equation  21,  Equation  22  becomes 


dfi 


ds  I  n( P)  \ 
dp  V  "  n(P')  / 


Now,  to  first  order  in  dp  we  have 

n(P')  =  n(P )  +  Vn(P)-Gdp. 
Substituting  into  Equation  23  gives 


,_(l  +  !-£L« tef1 

'  n(P)  / 


Again  to  first  order  in  dp  we  have 


Vn(P)  •  udp  \_l 
n(P)  ' 


Vn(P)  -{idp 
n(P) 


(23) 


(24) 
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Substituting  into  the  previous  expression  and  simplifying,  we  finally  obtain 


d/? 


Vn- u 

- ds, 


(25) 


where  Vn  and  n  are  henceforth  understood  to  be  evaluated  at  P. 

Equation  25  shows  that  the  rotation  of  It  towards  u  is  proportional  to  Vn  •  ft.  If  v  is  another 
unit  vector,  perpendicular  to  bothii  and  ft  (i .e..  v  =  kxu),  then  we  can  similarly  show  that  the 
rotation  of  ft  towards  v  is  proportional  to  Vn  •  v.  But  if  ft  is  perpendicular  to  both  It  and  ft,  then 
(Figure  1)  ft  must  also  be  perpendicular  to  Vn;  therefore,  Vn  •  v  =  0,  and  there  is  no  rotation  of  It 
toward  v.  In  other  words,  ft(Q)  in  Figure  1  indeed  lies  in  the  plane  offt(P)  and  Vn(P) 


Equation  25  is  the  principal  result  of  this  section.  This  equation  gives  the  infinitesimal 
angle  dp  through  which  the  unit  ray  vector  it  will  turn  in  the  next  infinitesimal  distance  ds  along 
the  ray.  In  interpreting  Equation  25,  we  must  remember  that  ft  lies  in  the  plane  of  £  and  Vn. 
making  a  right  angle  with  1c  and  an  acute  angle  with  Vn\  furthermore,  the  positive  sense  of  d/J  is 
from  ft  toward  ft.  The  unit  ray  vector  thus  turns  toward  regions  of  higher  n  in  the  plane  defined 
by  itself  and  Vn. 

The  local  turning  rate  of  the  ray  is  formally  defined  to  be  dp/ds,  and  is  immediately 
calculable  from  Equation  25.  Figure  1  shows  that  ds/d/J,  the  reciprocal  of  the  local  turning  rate,  is 
the  local  radius  of  curvature  of  the  ray. 


THE  RAY  EQUATIONS  FOR  n  =  n(|r|) 


We  shall  now  use  the  equation  for  the  turning  rate  of  the  unit  ray  vector.  Equation  25,  to 
derive  a  set  of  differential  equations  from  which  the  ray  trajectory  can  be  calculated  in  the  case 
that  the  index  of  refraction  n  at  any  point  r  is  a  function  only  of  r  This  case  serves  as  a  simplified 
model  of  Earth’s  atmosphere,  with  the  origin  at  the  center  of  the  Earth  and  r  bounded  below  bv 
the  radius  of  the  Earth. 


Let  P  be  any  point  in  the  medium  with  position  vector  r  =  rft  relative  to  the  origin  O,  and 
let  ft  be  the  unit  ray  vector  at  P  See  Figure  2.  Let  ft  be  any  constant  unit  vector  in  the  plane  of  ft 
and  ft,  and  let  6  be  the  angle  between  ft  and  ft  Finally,  let  ft  be  the  unit  vector  in  the  plane  of  ft 
and  ft  that  is  perpendicular  to  ft  and  in  the  direction  of  increasing  Q,  and  let  e  be  the  angle  between 
ft  and  ft.  We  may  regard  e  as  the  "elevation  angle”  of  ft.  The  "state”  of  the  ray  can  evidently  be 
defined  by  the  three  variables  r,  d  and  s.  Essentially,  r  and  6  tell  us  where  the  ray  is,  while  e  tells 
us  where  the  ray  is  going.  As  shown  in  Figure  3,  movement  along  the  ray  through  an 
infinitesimal  distance  ds  will  bring  us  to  a  point  Q  whose  position  vector  relative  to  O  is 
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Let  N  be  the  intersection  of  the  line  OQ  and  the  line  through  P  parallel  to  0.  Since  ds  is 
infinitesimally  small,  the  triangle  PQN  may  be  regarded  as  a  right  triangle  with  hypotenuse  PQ 
=  ds  and  sides  NQ  =  dr  and  PN  =  rdd.  Thus,  from  Figure  3  we  have 

sine  =  (dr)/(ds)  and  cose  =  (rd0)/(ds), 

whence 


dr  =  sineds,  (26) 

and 

d0  =  r-1cose  ds.  (27) 


Equations  26  and  27  give  the  respective  changes  in  r  and  0  as  we  move  an  infinitesimal 
distance  ds  along  the  ray  from  P.  We  have  only  to  calculate  the  concomitant  change  in  e.  For  this 
we  must  use  the  turning  rate  equation  (25). 

The  assumption  n(r)  =  n(r)  implies  that 

Vn  =  n'(r)ft.  (28) 

In  Figure  4  we  show  Vn  ,  along  with  the  unit  vector  &  that  lies  in  the  plane  of  £  and  Vn  making  a 
right  angle  with  k  and  an  acute  angle  with  Vn  .  We  also  show  in  Figure  4  the  angle  /?  that  ft 
makes  with  a  fixed  direction  in  space;  we  have  taken  that  fixed  direction  to  be  —ft  in  order  to 
comply  with  the  requirement  that  an  increase  in  P  correspond  to  a  rotation  of  ft  toward  u  (See 
discussion  following  Equation  25.)  Since  the  geometry  of  Figure  4  implies  that  the  angle  between 
ft  and  ft  is  equal  to  e,  we  have 

Vn-vi  =  n'(r)  r-u  =  n'{r)c ose. 

Therefore,  Equation  25  gives 

dp  =  [n'(r)/n(r)]  cose  ds.  (29) 

From  the  geometry  of  the  triangle  in  Figure  4  containing  the  two  angles  0  and  P,  it  is  seen 

that 

6  +  P  +  (nr/2  —  e)  =  n.  (30) 

Therefore, 

£  =  0  +  P  —  nl  2, 

whence 

de  =  d0  +  d  p.  (31) 

Substituting  on  the  right  from  Equations  27  and  29,  we  deduce 

de  =  (r_1  +  n'(r)/n(r)]coseds.  (32) 

Equation  32  is  the  desired  equation  relating  de  to  ds. 
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In  summary,  if  current  values  of  the  state  variables  r,  0  and  e  of  the  ray  are  known,  then 
displacement  along  the  ray  by  an  infinitesimal  distance  ds  will  change  these  state  variables  by 
the  respective  infinitesimal  amounts  [cf.  Equations  26, 27,  and  32] 

(dr  =  sine  ds,  (33a) 

d0  =  r-1coseds,  (33b) 


|de  =  [r  1  +  n'(r)/n(r)]coseds.  (33c) 

Equations  33  constitute  the  main  result  of  our  analysis.  They  imply  that  r,  0,  and  e  can  be 
regarded  as  functions  of  the  ray  length  s,  which  functions  may  be  found  by  solving  the  set  of 
coupled,  first-order,  ordinary  differential  equations 

dr/ds  =  sine,  (34a) 

<  d0/ds  =  r-1  cose,  (34b) 

de/ds=  [r-1  +  ra'(r)/n(r)]  cose,  (34c) 

subject  to  the  specified  initial  conditions  r  =  rQl  0  =  0Q  and  e  =  e0  at  s  =  0. 

Actually,  we  can  write  down  three  more  equivalent  sets  of  differential  equations  by 
successively  taking  r,  0  and  e  as  the  independent  variable  instead  of  s.  For  example,  if  we  solve 
Equations  33  for  dr,  ds,  and  de  in  terms  of  d0,  we  may  deduce  the  set  of  differential  equations 

dr/d0  =  r~ 1  tane,  (35a) 

.  ds/d0  =  rt cose,  (35b) 

de/d0  =  1  +  r-1n'(r)/n(r).  (35c) 

The  advantage  of  using  s  as  the  independent  variable  is  that  the  resulting  differential  equations 
34  are  well  behaved  for  r  bounded  away  from  zero.  By  contrast,  for  example,  Equations  35a  and 
35b  be  <e  unbounded  as  ±  n/2. 


SOLVING  THE  RAY  EQUATIONS 

Keeping  in  mind  that  the  ray  equations  33  and  their  associated  o  d  e.  forms  are 
approximately  valid  only  for  light  rays  of  sufficiently  high  frequency,  let  us  now  consider  how  to  go 
about  solving  these  equations. 


Case  n(r)  =  Constant 

If  n(r)  is  independent  of  r,  then  the  second  term  in  brackets  in  Equation  33c  vanishes,  and 
we  have 
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where  the  second  equality  has  invoked  Equation  33b.  Integrating,  we  get 


e  -  Eo  =  9  -  d0- 


(37) 


Now  dividing  Equation  33a  by  Equation  33b,  and  then  substituting  for  £  from  Equation  37,  we  get 


dr  _  rsin(d-d0  +  sQ) 
dd  cosld  —  9Q  +  £^) 


dr 

r 

Equation  38  is  easily  integrated  to  give 


d  cos(0-0o+£o) 
cos(0-0o  +  £o) 


(38) 


whence 


/  COS (d-d0  +  E0)  '  /  COS£0 

\  COS(9Q-9Q  +  £^  /  \  COS(0-0Q-f  £q) 


rcos(0-0o  +  £o)  =  r0cosco.  (39) 

Equation  39  is  just  the  polar  equation  of  a  straight  line,  which  is  precisely  the  ray  equation  we 
expect  when  the  index  of  refraction  is  constant. 


Case  n(r)  =  A/r. 

If  n(r )  =  Air,  where  A  is  a  constant,  then  we  have  n'(r)  =  —  A/r2.  so  that 

n'(r)/n(r)  =  (-Air2)  I  (A/r)  =  -r~l.  (40) 

In  that  case,  the  right  side  of  Equation  33c  vanishes  everywhere,  so  that  Equation  33c  integrates 
to 


£  =  £0'  <41> 

Thus,  the  trajectory  is  such  that  the  elevation  angle  of  the  unit  ray  vector  stays  constant.  If  £0  = 
0,  then  Equations  34a  and  34b  are  easily  integrated  to  give  r  =  rQ  and  6  =  (9Q  4-  s/rQ,  which  is 
the  equation  of  a  circle.  If  £Q  *  0,  then  integration  of  Equations  34a  and  34b  give 

r  =  r0  +  ssin£0,  (42a) 

6  =  6q  +  cot£0  log[l  +(s/r0)sin£0|.  (42b) 
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These  equations  describe  a  curve  that  spirals  either  inward  or  outward,  according  to  whether  eQ  is 
negative  or  positive. 


Numerical  Integrations 

For  even  moderately  realistic  functions  n(r),  usually  we  must  proceed  with  a  computer- 
oriented  numerical  calculation.  The  simplest,  though  not  the  most  efficient,  method  is  to 
recursively  apply  Equations  33  with  the  differentials  dr,  dd,  de,  and  ds  replaced  by  "small  but 
finite  increments”  A r,  Ad,  Ae,  and  As,  respectively.  Figure  5  summarizes  this  method  of 
constructing  the  ray;  essentially,  it  is  the  Euler  method  of  integrating  the  differential  equations 
34,  subject  to  the  initial  conditions  r  =  rQ,  d  =  dQ,  £  =  eQ  at  s  =  0.  The  key  here  is  to  choose  As  so 
small  that  reducing  it  further  does  not  substantially  alter  the  results.  Of  course,  the  smaller  As 
is,  the  longer  it  will  take  to  trace  the  ray  through  a  given  length  s.  The  stopping  condition  in  step 
5  of  Figure  5  will  generally  be  that  one  of  the  variables  r,  d,  £,  or  s  has  reached  some 
predetermined  value;  or,  it  might  be  that  a  predetermined  value  has  been  reached  by  the  travel 
time  t  along  the  ray,  the  increment  of  which  is  given  by 

At  =  (n(r)/c)  As,  (43) 

since  dn(r)  is  the  wave  velocity  at  r. 

A  more  efficient  way  of  solving  the  set  of  differential  equations  (34)  would  be  to  simply  use 
a  packaged  computer  code.  An  example  is  the  IMSL  subroutine  DVERK,  which  is  available  at 
most  computer  centers.  Subroutine  DVERK  uses  a  high-order  Runge-Kutta  method,  which  yields 
a  much  more  accurate  estimate  of  the  solution  functions  Ks),  dis),  and  e(s)  for  a  given  step  size  As, 
as  compared  to  the  simpler  Euler  method  of  Figure  5. 

Regardless  of  what  numerical  method  is  selected  to  integrate  Equations  34,  a  key  user 
input  is  the  index  of  refraction  function,  n(r).  At  each  As  step  the  logarithmic  derivative  n’(r)/n(r) 
must  be  computed  anew,  either  from  some  assumed  formula  or  else  from  tabulated  data. 

A  strong  consistency  check  on  any  numerical  integration  of  Equations  34  can  be  had  by 
periodically  testing  how  well  the  conservation  law  in  Equation  44  is  satisfied. 

A  Conservation  Law 

We  shall  now  show  that  Equations  33  imply  that 

r  n(r)  cose  =  rQ  n(r0)  coseQ.  (44)  ^ 

In  other  words,  the  quantity  on  the  left  side  of  Equation  44  is  constant  along  any  ray  trajectory. 

Dividing  Equation  33c  by  33a,  wc  obtain 

de  [r-1  +  n'(r)/ n(r)l  cose 
dr  sine 
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Rearranging,  we  get 


or  equivalently, 


Integrating  gives 


sinrdc  _  dr  n'(r)dr 

cose  r  n(r) 


d( cose) 
cose 


dr  d(rt(r)> 

r  n(r) 


This  last  equation  can  also  be  written 


,  f  r  n(r)  cose  \ 

In  - )  =  0 

\  r  n(r  )  cose„  / 


ln(l), 


from  which  Equation  44  immediately  follows. 

Equation  44  can  be  viewed  as  a  sort  of  generalization  of  Snell’s  law.  The  latter  states  that  if 
a  light  ray  crosses  a  planar  surface  at  which  the  index  of  refraction  changes  discontinuously  from 
nj  to  n^,  then  the  angles  et  and  e2  between  the  surface  and  the  rays  on  either  side  satisfy  the 
equation  rijcoscj  =  n2 cose2.  Although  this  last  equation  can  in  fact  be  derived  from  Equation  44 
by  using  a  limiting  argument,  one  cannot  legitimately  proceed  the  other  way  around  and  derive 
the  conservation  law.  Equation  44,  or  the  ray  equations  (33),  from  Snell’s  law  instead  of  from 
Maxwell’s  equations.  The  reason  for  this  is  that  Snell’s  law  applies  to  plane  waves  of  any 
frequency  at  a  planar  discontinuity  in  the  index  of  refraction,  whereas  the  conservation  law, 
Equation  44  and  ray  equations  (33)  apply  to  distorted  plane  waves  of  high  frequency  only  in  a 
medium  where  the  index  of  refraction  is  a  continuous  (indeed,  differentiable)  function  of  r. 


"Flat  Earth”  Formulas 

It  often  happens  that  we  are  interested  in  ray  trajectories  in  the  Earth’s  atmosphere  over 
distances  that  are  extremely  small  compared  to  the  Earth’s  radius  RE.  In  that  case  it  is  often 
permissible  to  assume  that  the  ray  is  propagating  in  an  xz  plane,  with  the  z  axis  pointing  ”up”  and 
the  x  axis  "horizontal”  (  and  in  the  plane  of  the  z  axis  and  the  initial  unit  ray  vector).  The  r-only 
dependence  of  the  index  of  refraction  n  easily  translates  into  a  z-only  dependence.  As  shown  in 
Figure  6,  for  sufficiently  small  6  we  can  define  x  and  z  in  terms  of  r  and  9  by 

z  »  r  -  Re,  and  x  »  rsinfl. 
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It  follows  that 


d z  ~  dr, 


and 


dx  =  dr  sin0  +  r  cos©  d0, 

or,  since  6  =»  0  for  the  relatively  small  propagation  distances  of  interest. 


(45a) 


dx  »  rd0.  (45b) 

Because  of  Equations  45,  we  may  replace  dr  in  Equation  33a  with  dz,  and  rdO  in  Equation 
33b  with  dx.  Furthermore,  we  may  assume  that  RE  is  so  large  that  r~ 1  =  ( RE  +  z )  ~ 1  is  negligible 
compared  to  n'(r)/n(r)  =  n'(z)/n(z )  in  Equation  33c.  We  thus  conclude  that  the  equations 
describing  light  ray  trajectories  in  a  "flat  Earth”  atmosphere,  in  which  the  index  of  refraction  n  is 
a  function  only  of  the  altitude  z,  are: 


dz  =  sine  ds, 

(46a) 

dx  =  cose  ds. 

(46b) 

de  =  [n'(2)/n(z)]  cose  ds. 

(46c) 

Notice  that  Equations  46a,b  imply  that  dz/dx  =  tane,  showing  that  the  elevation  angle  £  is  just 
the  conventional  "slope  angle”  of  the  ray  trajectory  in  the  xz  plane. 

As  with  Equations  33,  one  can  proceed  to  write  down  from  Equations  46  four  different  sets 
of  three  coupled  differential  equations,  one  set  for  each  choice  of  z,  x,  e,  or  s  as  the  independent 
variable.  And  also  as  before,  the  set  with  s  as  the  independent  variable  is  usually  the  best 
behaved. 


The  conservation  law  in  Equation  44  now  reads 

( 2  +  Re*  n(z )  cose  *»  constant. 

But  since  RE  >  z,  we  can  drop  the  z  from  the  lead  factor  and  then  absorb  RE  into  the  constant  on 
the  right.  Thus,  the  "flat  Earth”  conservation  law  reads 

n(z)  cose  =  n(z0)  cos£Q.  (47) 

Indeed,  it  is  easy  to  derive  Equation  47  directly  from  Equations  46:  We  merely  divide  Equation 
46c  by  Equation  46a,  separate  variables,  and  then  integrate,  just  as  we  did  in  deriving  Equation 
44  from  Equations  33.  Equation  47  can  evidently  be  used  as  a  consistency  check  on  any  numerical 
integration  of  Equations  46. 


FIGURE  1.  Geometry  for  deriving  from  the  eikonal  equation  an 
expression  for  the  infinitesimal  angle  d/?  through  which  the  unit  ray  vector 
^  rotates  in  moving  an  infinitesimal  distance  ds  along  the  ray. 
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FIGURE  2.  Geometry  showing  the  variables  r,  d  and  e  that  define  the 
instantaneous  state  of  a  ray  passing  through  a  point  P  in  the  direction  ft. 
The  origin  of  s,  the  path  length  along  the  ray,  is  arbitrary. 
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1.  Initialize:  Set  r  =  rQ,  0  =  6^ ,  e  =  s0  and  s  =  0. 

I 

2.  Choose  a  small,  positive  value  for  As. 


J 

r3-  Calculate,  for  the  current  values  of  r  and  e,  the  quantities: 


A  r  =  sine  •  As, 


A 0  =  r  1  cose  •  As, 


Ae  =  [r  1  +  n'(r)/n(r)]  cose  ■  As. 


4.  Update: 


r  =  r  +  Ar, 
0  =  0  +  Ad, 
e  =  e  +  Ae, 
s  =  s  +  As. 


L5.  Record  the  current  values  of  r,  0,  e,  ands.  Stop  if  the 
termination  condition  is  satisfied;  otherwise,  go  to  Step  3. 


FIGURE  5.  The  simple  Euler  method  for  numerically  integrating 
Equations  34,  given  that  at  s  =  0  we  have  r  =  rQt  0=0oand  e  =  eQ  In  Step  2, 
As  should  be  chosen  small  enough  that  a  repeat  run  with  As  chosen  to  be 
only  half  as  large  should  not  significntly  alter  the  terminal  values  of  the 
variables. 
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FIGURE  6.  Geometry  for  deducing  the  form  of  the  ray  equations  in 
Earth’s  atmosphere  in  situations  where  the  curvature  of  the  Earth’s  surface 
is  negligible. 


